!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! 
!! Turbulence closures
!!
!! \n
!!
!! This module provides a collection of turbulence closures for the
!! (weakly) compressible Navier-Stokes equations. This module is
!! complicated because it has to consider the following aspects.
!! <ul>
!!  <li> Each turbulent model requires specific pointwise values, both
!!  for variables and their gradients. While the pointwise values for
!!  the variables can be computed easily from the main prognostic and
!!  diagnostic quantities used in \c mod_dgcomp_rhs, each gradient
!!  requires a global computation and must be computed in \c
!!  mod_dgcomp_rhs. Ideally, depending on the chosen turbulence model,
!!  only the relevant gradients should be computed.
!!  <li> Some turbulent models require element-averaged quantities,
!!  which are better computed in a specific element loop. Moreover,
!!  such quantities should be precomputed and stored for all the
!!  elements in the mesh for at least two reasons:
!!  <ul>
!!   <li> during the side loop, the element averaged quantities from
!!   the two neighbouring elements are required;
!!   <li> one could think to update such quantities only at selected
!!   time levels, for efficiency reasons.
!!  </ul>
!!  <li> The element averaged quantities typically involve some
!!  gradients, thereby introducing a coupling between the two previous
!!  requirements.
!!  <li> It is often desirable to obtain detailed diagnostics which
!!  are specific to a given turbulence model; this implies that this
!!  module must interact with \c mod_dgcomp_statistics.
!!  <li> The turbulence model must be evaluated for a generic state \c
!!  uuu, without any assumption about previous evaluations.
!!  <li> The turbulence model might require additional prognostic
!!  variables, such as the turbulent kinetic energy related variables
!!  in the \f$\kappa-\epsilon\f$ and \f$\kappa-\omega\f$ models.
!! </ul>
!! To meet all these requirements, we use here a three-layer,
!! "sandwich" layout: a base module \c mod_turb_flux_base, defining
!! the general layout, a collection of modules specific to various
!! turbulence models, and the present interface module. The turbulence
!! module is described by three types: the module itself, the module
!! diagnostic variables (typically, gradients and coefficients), and
!! the module prognostic variables. This distinction is based on the
!! general layout adopted for ODE problems in \c mod_time_integrators.
!!
!! The typical work flow is
!! \code{.f90}
!!  call turb_mod%init(vars) ! allocate the working arrays
!!
!!  do it=1,nt ! time loop
!!
!!    ! Compute the gradients: implemented in mod_dgcomp_rhs with calls
!!    ! to turb_mod%compute_grad_diags
!!
!!    ! Compute the element averaged constants, as well as other
!!    ! diagnostic quantities
!!    call turb_mod%compute_coeff_diags(vars)
!!
!!    ! Compute the turbulent fluxes: on the fly for elements and sides
!!    ! If required, update the specific statistics passing the
!!    ! optional argument sc
!!    do ie=1,ne ! element loop
!!      call turb_mod%flux(ie,vars,sc)
!!    enddo
!!    do is=1,ns ! side loop
!!      do ie=1,2 ! loop on the two connected elements
!!        call turb_mod%flux(ie,vars,sc)
!!      enddo
!!    enddo
!!
!!  enddo
!! \endcode
!!
!! The turbulent fluxes are collected in a three dimensional array
!! \f${\bf F}_{turb}\f$ as follows:
!! \f{displaymath}{
!!  \begin{array}{ccl}
!!  {\bf F}_{turb}(\,:\,,\phantom{2:}1\phantom{+d},l) & = &
!!  \underline{\mathcal{F}}^d_E(x_l)
!!  \\[1mm]
!!  {\bf F}_{turb}(\,:\,,2:1+d,l) & = &
!!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}}(x_l)
!!  \end{array}
!! \f}
!! where \f$x_l\f$, \f$l=1,\ldots,N_{points}\f$ are the coordinates of
!! the points where the flux is required.
!!
!! \note The functions provide also some optional arguments which can
!! be used to recover intermediate values used in the computation of
!! the turbulent fluxes.
!!
!! \section uniform_expansion Uniform expansion
!!
!! It is useful to derive the relation between the symmetric velocity
!! gradient and the velocity divergence for the case of a uniform
!! expansion, since sometimes the dissipative fluxes are defined so
!! that they are zero for such a flow. Let us denote by \f$d\f$ the
!! spatial dimension and let us consider a flow characterized by:
!! \f{displaymath}{
!!  \nabla\underline{u} = \eta\,\mathcal{I}
!! \f}
!! for some expansion velocity \f$\eta\f$. We have that
!! \f{displaymath}{
!!  \nabla\underline{u}+\nabla\underline{u}^T = 2\eta\, \mathcal{I} =
!!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}.
!! \f}
!! Hence, to guarantee that a uniform expansion does not generate any
!! dissipative stress, we can express the dissipative fluxes as a
!! function of the deviatoric deformation velocity
!! \f{displaymath}{
!!  \mathcal{S}_D = \nabla\underline{u}+\nabla\underline{u}^T -
!!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}.
!! \f}
!! 
!! \subsection Example
!!
!! Consider the flow
!! \f{displaymath}{
!!  u_1 = \eta x_1,\qquad u_2 = \eta x_2, \qquad u_3 = 0.
!! \f}
!! If this flow is seen as a two-dimensional flow, \f$\mathcal{S}_D\f$
!! vanishes, however if it is regarded as a three-dimensional flow we
!! have
!! \f{displaymath}{
!!  \mathcal{S}_D = \frac{\eta}{3}\left[ 
!!   \begin{array}{ccc}
!!     2 & 0 & 0 \\
!!     0 & 2 & 0 \\
!!     0 & 0 & -4 \\
!!   \end{array}
!!   \right].
!! \f}
!<----------------------------------------------------------------------
module mod_turb_flux

!-----------------------------------------------------------------------

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_turb_flux_base, only: &
   mod_turb_flux_base_constructor, &
   mod_turb_flux_base_destructor,  &
   mod_turb_flux_base_initialized, &
   c_turbmod, c_turbmod_progs, c_turbmod_diags, &
   t_turb_diags
 
 use mod_viscous_flux, only: &
   mod_viscous_flux_constructor, &
   mod_viscous_flux_destructor,  &
   mod_viscous_flux_initialized, &
   t_viscous_flux

 use mod_smagorinsky_flux, only: &
   mod_smagorinsky_flux_constructor, &
   mod_smagorinsky_flux_destructor,  &
   mod_smagorinsky_flux_initialized, &
   t_smagorinsky_flux

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid 

 use mod_dgcomp_testcases, only: &
   mod_dgcomp_testcases_initialized, &
   t_phc, phc, coeff_visc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_turb_flux_constructor, &
   mod_turb_flux_destructor,  &
   mod_turb_flux_initialized, &
   c_turbmod, c_turbmod_progs, c_turbmod_diags, t_turb_diags, &
   t_viscous_flux, t_smagorinsky_flux

 private

!-----------------------------------------------------------------------

! Module types and parameters

! Module variables

 ! public members
 logical, protected ::               &
   mod_turb_flux_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_turb_flux'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_turb_flux_constructor(base, grid, k_dyn)
  type(t_base), intent(in) :: base
  type(t_grid), intent(in) :: grid
  integer, intent(in)      :: k_dyn
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
(mod_dgcomp_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_turb_flux_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   call mod_turb_flux_base_constructor()
   call mod_viscous_flux_constructor()
   call mod_smagorinsky_flux_constructor()

   mod_turb_flux_initialized = .true.
 end subroutine mod_turb_flux_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_turb_flux_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_turb_flux_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call mod_smagorinsky_flux_destructor()
   call mod_viscous_flux_destructor()
   call mod_turb_flux_base_destructor()
   mod_turb_flux_initialized = .false.
 end subroutine mod_turb_flux_destructor

!-----------------------------------------------------------------------


!!-----------------------------------------------------------------------
!
! !> Smagorinsky-Richardson isotropic turbulent viscosity
! !!
! !! In this subroutine we consider the model discussed in <a
! !! href="http://onlinelibrary.wiley.com/doi/10.1111/j.2153-3490.1962.tb00128.x/pdf">
! !! [Lilly, 1962]</a> and <a
! !! href="http://journals.ametsoc.org/doi/pdf/10.1175/1520-0493%281983%29111%3C2341%3AACMFTS%3E2.0.CO%3B2">
! !! [Durran, Klemp, 1983]</a>. In practice, this amounts to a
! !! Smagorinsky eddy viscosity corrected according to the local
! !! Richardson number.
! !!
! !! We set
! !! \f{displaymath}{
! !!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}} =
! !!  -\rho\nu_U\mathcal{S}, \qquad
! !!  \underline{\mathcal{F}}^d_E = -c_p\rho\nu_E \nabla T.
! !! \f}
! !! with
! !! \f{displaymath}{
! !!  \mathcal{S} = \nabla\underline{u} + \nabla\underline{u}^T -
! !!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}.
! !! \f}
! !! The turbulent viscosities are
! !! \f{displaymath}{
! !!  \nu_U = k^2 \Delta x\,\Delta z\, \frac{|\mathcal{S}|}{\sqrt{2}}
! !!  \sqrt{\max\left(1-\frac{Ri}{Pr}\,,\sigma_{min}\right)}, \qquad
! !!  \nu_E = Pr^{-1}\nu_U,
! !! \f}
! !! where \f$k=0.21\f$, \f$Pr=\frac{1}{3}\f$ and
! !! \f{displaymath}{
! !!  |\mathcal{S}| = \sqrt{\sum_{i,j}\mathcal{S}^2_{ij}} , \qquad Ri =
! !!  2 \frac{\frac{g}{\theta} \frac{\partial\theta}{\partial
! !!  z}}{|S|^2}.
! !! \f}
! !! The coefficient \f$\sigma_{min}\f$ should be zero, which would
! !! imply the possibility of completely suppressing turbulent mixing.
! !! A value larger than zero is often useful for numerical stability.
! pure subroutine smagrich_flux(fem , d,h,x,rho,p,uu , ti , turb_diags)
!  integer, intent(in) :: d
!  real(wp), intent(in) :: h
!  real(wp), intent(in) :: x(:,:)
!  real(wp), intent(in) :: rho(:), p(:)
!  real(wp), intent(in) :: uu(:,:)
!  type(t_turb_input), intent(in) :: ti
!  real(wp), intent(out) :: fem(:,:,:)
!  type(tttt_turb_diags), intent(out), optional :: turb_diags
!
!  logical, parameter :: local_ri = .false.
!
!  real(wp), parameter :: &
!    k  = 0.21_wp,       &
!    pr = 1.0_wp/3.0_wp, &
!    smin = 0.001_wp
!  integer :: l, id
!  real(wp) :: theta(size(x,2)), ldivu, sym(d,d), ri, mu
!
!   ! Compiler bug: this should not be necessary, however see
!   ! http://software.intel.com/en-us/forums/showthread.php?t=82775&o=a&s=lr
!   ! and bug DPD200169598 for ifort.
!   if(present(turb_diags)) then
!     if(allocated(turb_diags%s_abs)) deallocate(turb_diags%s_abs)
!     if(allocated(turb_diags%ri   )) deallocate(turb_diags%ri   )
!     if(allocated(turb_diags%nu   )) deallocate(turb_diags%nu   )
!     if(allocated(turb_diags%diss )) deallocate(turb_diags%diss )
!     if(allocated(turb_diags%s2_e )) deallocate(turb_diags%s2_e )
!   endif
!   ! end of compiler bug
!
!   if(present(turb_diags)) then
!     allocate(turb_diags%s_abs(size(x,2)))
!     allocate(turb_diags%ri   (size(x,2)))
!     allocate(turb_diags%nu   (size(x,2)))
!     allocate(turb_diags%diss (size(x,2)))
!     allocate(turb_diags%s2_e (size(x,2)))
!   endif
!
! !!
! !! The potential temperature is computed as
! !! \f{displaymath}{
! !!  \theta = \left(\frac{p_s}{R}\right)^\kappa \frac{T}{(\rho
! !!  T)^\kappa}.
! !! \f}
! 
!   theta = p/(rho*phc%rgas) * (phc%p_s/p)**phc%kappa
!   do l=1,size(x,2)
!     ! momentum flux
!     ldivu = (-2.0_wp/real(d,wp)) * 0.5_wp*tr(ti%gsuu(:,:,l))
!     sym = ti%gsuu(:,:,l) ! symmetric gradient
!     do id=1,d
!        sym(id,id) = sym(id,id) + ldivu
!     enddo
!
!     if(local_ri) then
!       ! Richardson number multiplied by |sym|^2 to avoid 0 division
!       ri = 2.0_wp * phc%gravity*ti%gth(d,l)/theta(l)
!       mu = rho(l) * (k**2/sqrt(2.0_wp)) * h**2  &
!             * sqrt(max( sum(sym**2)-ri/pr , 0.0_wp ))
!     else
!       ! all the terms including derivatives are averaged over the
!       ! gradient to compensate the numerical instability of 
!       ri = 2.0_wp * phc%gravity*ti%gzth_e/theta(l)
!       mu = rho(l) * (k**2/sqrt(2.0_wp)) * h**2  &
!             * sqrt(max( ti%s2_e-ri/pr , smin*ti%s2_e ))
!     endif
!
!     fem(:,2:1+d,l) = -mu * sym
!
!     ! energy flux
!     fem(:,  1  ,l) = -phc%cp * mu/pr * ti%gtt(:,l)
!
!     ! tracers
!     fem(:,2+d: ,l) = -mu/pr * ti%gcc(:,:,l)
!
!     if(present(turb_diags)) then
!       turb_diags%s_abs(l) = sqrt(sum(sym**2))
!       turb_diags%ri   (l) = ri ! more precisely:  |S|^2 * ri = 2*N^2
!       turb_diags%nu   (l) = mu/rho(l)
!       turb_diags%diss (l) = -sum(fem(:,2:1+d,l)*sym) ! use sym^T = sym
!       turb_diags%s2_e (l) = sqrt(ti%s2_e)
!     endif
!   enddo
!
! end subroutine smagrich_flux
!
!!-----------------------------------------------------------------------
!
! !> Anisotropic model.
!!  !!
!!  !! In this subroutine we follow the general framework of <a
!!  !! href="http://www.sciencedirect.com/science?_ob=MImg&_imagekey=B6TYJ-4B82P8M-3M-1&_cdi=5620&_user=2620285&_pii=S0898122103900149&_origin=gateway&_coverDate=08%2F31%2F2003&_sk=999539995&view=c&wchp=dGLzVlz-zSkWb&md5=bd04c1f14bad52d9b4db33591284d7ee&ie=/sdarticle.pdf">
!!  !! [Abb&agrave;, Cercignani, Valdettaro, 2003]</a> to prescribe a
!!  !! turbulent diffusion in the vertical direction as discussed in <a
!!  !! href="http://onlinelibrary.wiley.com/doi/10.1111/j.2153-3490.1962.tb00128.x/pdf">
!!  !! [Lilly, 1962]</a>. In practice, this amounts to a Smagorinsky eddy
!!  !! viscosity corrected according to the Richardson number and
!!  !! restricted to the vertical direction.
!!  !!
!!  !! Let us first consider the momentum flux and let us write, for
!!  !! simplicity, \f$\tau_{ij} = \left(
!!  !! \underline{\underline{\mathcal{F}}}^d_{\underline{U}}
!!  !! \right)_{ij}\f$. The idea is to write
!!  !! \f{displaymath}{
!!  !!  \tau_{ij} = -2 \tilde{\Delta}^2|S|B_{ijrs} S_{rs},
!!  !! \f}
!!  !! where \f$S = \nabla\underline{u} + \nabla\underline{u}^T\f$ and we
!!  !! have kept the isotropic component. The diffusion tensor is now
!!  !! assumed to have the following form:
!!  !! \f{displaymath}{
!!  !!  B_{ijrs} = \sum_{\alpha,\beta} C_{\alpha\beta}\, a_{i\alpha}\,
!!  !!  a_{j\beta}\, a_{r\alpha}\, a_{s\beta}
!!  !! \f}
!!  !! where \f$a\f$ denotes the normal matrix whose columns are \f$d\f$
!!  !! orthonormal vectors \f$\underline{a}_\alpha\f$ and
!!  !! \f$C_{\alpha\beta}\f$ is a symmetric tensor. If we now consider a
!!  !! reference system aligned with the unit vectors
!!  !! \f$\underline{a}_\alpha\f$, so that \f$a_{ij}=\delta_{ij}\f$, we
!!  !! obtain
!!  !! \f{displaymath}{
!!  !!  \tau_{ij} = -2 \tilde{\Delta}^2|S|C_{ij} S_{ij} \qquad {\rm
!!  !!  (no\,\,sum)}.
!!  !! \f}
!!  !! To determine the values \f$C_{ij}\f$, let us consider the
!!  !! hypothesis that the momentum diffusion only takes place along the
!!  !! vertical direction, which we denote by \f$z\f$, i.e.
!!  !! \f{displaymath}{
!!  !!  \tau_{\mu\nu} = \delta_{z\}
!!  !! \f}
! !<
!! pure subroutine aniso_flux(fem , d,x,rho,uu,guu,gtt)
!!  integer, intent(in) :: d
!!  real(wp), intent(in) :: x(:,:)
!!  real(wp), intent(in) :: rho(:)
!!  real(wp), intent(in) :: uu(:,:), guu(:,:,:)
!!  real(wp), intent(in) :: gtt(:,:)
!!  real(wp), intent(out) :: fem(:,:,:)
!!
!!  integer :: l, id
!!  real(wp) :: nu(2,size(x,2)), ldivu, sym(d,d)
!! 
!!   ! problem coefficients
!!   nu = coeff_visc( x , uu )
!!
!!   do l=1,size(x,2)
!!     ! momentum flux
!!     ldivu = (-2.0_wp/3.0_wp) * tr(guu(:,:,l)) ! lambda*div(u)
!!     sym = guu(:,:,l)+transpose(guu(:,:,l)) ! symmetric gradient
!!     do id=1,d
!!       sym(id,id) = sym(id,id) + ldivu
!!     enddo
!!     fem(:,2:,l) = -rho(l)*nu(2,l)*sym
!!
!!     ! energy flux
!!     fem(:,1 ,l) = -rho(l)*nu(1,l)*phc%cp*gtt(:,l) &
!!                   + matmul( uu(:,l) , fem(:,2:,l) )
!!   enddo
!!
!! contains
!!
!!  pure function tr(a)
!!   real(wp), intent(in) :: a(:,:)
!!   real(wp) :: tr
!!   integer :: i
!!
!!    tr = a(1,1)
!!    do i=2,size(a,1)
!!      tr = tr + a(i,i)
!!    enddo
!!  end function tr
!!
!! end subroutine aniso_flux
!
!!-----------------------------------------------------------------------
!
! pure function tr(a)
!  real(wp), intent(in) :: a(:,:)
!  real(wp) :: tr
!  integer :: i
!
!   tr = a(1,1)
!   do i=2,minval(shape(a))
!     tr = tr + a(i,i)
!   enddo
! end function tr
!
!!-----------------------------------------------------------------------

end module mod_turb_flux

